Numerically Exact Long Time Behavior of Nonequilibrium Quantum Impurity Models 
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A Monte Carlo sampling of diagrammatic corrections to the non-crossing approximation is shown 
to provide numerically exact estimates of the long-time dynamics and steady state properties of 
nonequilibrium quantum impurity models. This 'bold' expansion converges uniformly in time and 
significantly ameliorates the sign problem that has heretofore limited the power of real-time Monte 
Carlo approaches to strongly interacting real-time quantum problems. The new approach enables the 
study of previously intractable problems ranging from generic long time nonequilibrium transport 
characteristics in systems with large onsite repulsion to the direct description of spectral functions 
on the real frequency axis in Dynamical Mean Field Theory. 



Numerical evaluation of diagrammatic perturbation 
series has played an important role in fields includ- 



statistical mechanics 
Recent developments 



ing quantum electrodynamics, |1| 
and condensed matter physics. [3 
have established diagrammatic Monte Carlo methods [8| 
as particularly powerful tools for the study of finite size 
clusters of interacting fermions coupled to a nonintcr- 
acting bath. These 'quantum impurity models' are used 
to model the physics of nanosystems coupled to leads 
Q, adsorption of atoms on surfaces [Io| and magnetic 
impurities in metals. 11 and are important as auxiliary 
problems in dynamical mean field theory approaches to 
infinite lattice correlated systems. [l2, 13 

Diagrammatic Monte Carlo methods have been very 
successful for equilibrium problems, where the Hamilto- 
nian H is partitioned as H = Hq + Hi the partition 
function Z at temperature T is expressed in an inter- 
action representation as Z = TrT T cxp — f^ T -Hi(r) 
and the exponential is expanded. The diagrammatic ex- 
pansion order needed to obtain reliable results grows as 
(Hi)/T, setting a lower limit on the temperatures which 
can be studied with fixed computational resources, but in 
practice the method works down to temperatures which 
are very low relative to the basic scales 14|. The use 



of Wick's theorem, which organizes fermion contractions 
into determinants, dramatically reducing the number of 
diagrams which must be sampled and eliminating one 
source of fermion sign problem, is crucial to the success 
of the procedure. 

In nonequilibrium problems the role of the partition 
function is played by the time evolution operator K ~ 
cxp [iffit] and the role of inverse temperature is played 
by the time interval t to be studied. The factor of i 
means that a straightforward expansion suffers from a 
phase problem, which in practice severely limits the dia- 
grammatic order which can be sampled and therefore the 
time intervals which can be studied. To date only rela- 
tively short times (up to ~ 2 — 3 times the hybridization 
scale) have been accessed [TBI - H^ ]. 

In analytic many-body theory, partial resummation 



techniques are often used to sum up specific classes of 
diagrams. If the resummation captures enough of the 
physics, one may hope that an expansion around it will 
converge rapidly. Motivatived by this idea, we formu- 
lated such an expansion for equilibrium properties of 
quantum impurity models (l9| . We found that while the 
method lead to a marked decrease in the mean pertur- 
bation order, the gain was in most cases offset by a sign 
problem, which is generic to this class of methods and 
arises from the absence of Wick's theorem for expansions 
about partially resummcd theories. 

In this paper we present a new formulation of a numer- 
ical expansion around a partial resummation, applicable 
to nonequilibrium quantum impurity models and there- 
fore also to time dependent dynamical mean field theory. 
We refer to this expansion as a 'bold' expansion but note 
that it differs from algorithms (also termed 'bold expan- 
sions') which sample higher-order self-energy diagrams 
based on numerically computed lower-order diagrams but 
employ no partial resummations [20I 21 1. Our method is 



based on a stochastic sampling of the full configuration 
space and is numerically exact. New aspects of our work 
include a treatment of the vertex corrections essential 
for the evaluation of expectation values, and most im- 
portantly a demonstration that in contrast to the bare 
expansion, the convergence of the bold expansion is uni- 
form in time: even the long time behavior is adequately 
characterized by finite orders of bold perturbation theory. 
We present results for the steady state density matrix, 
charge and magnetic relaxation times and current of a 
model quantum dot which demonstrate that the method 
allows access to unprecedentedly long times. 

We demonstrate the power of the method on the An- 
derson model with Hamiltonian 



H 



(e d + Ha) d\d a + Unfm + H hyh + H load (1) 



which describes a quantum dot with a single spin- 
degenerate orbital with correlation energy U hybridized 
to two leads labeled by a = L, R. The Hubert space of the 
impurity consists of four states: |0), | t), 4) and | t|). 
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H describes a magnetic field directed parallel to the spin 
quantization axis, H hyh = J2kaa \ykaod\c kaa + H.c] 
parametrizes the hybridization between the level and the 
leads and H\ ead describes the dynamics of the leads. Lead 
a is assumed to be in equilibrium at chemical poten- 
tial p a and temperature T a ; the presence of two leads 
allows for departures from equilibrium parametrized by 
Mi 7^ MJJ or Tl =/= Tr. An important parameter is the 
level width T = J2ka^ka^( £ka ~ Ma)- For our specific 
calculations we use the parametrization of Ref. [181 ] with 
Ml = "Mi? = ^7 2 > v = 10 = w c . 

We wish to compute time dependent expectation val- 
ues of operators O such as the dot charge (n) and spin 
(to) densities n-j- ± n± or the current flowing into the dot 
from (say) the left lead J L = i Y. ka [VkLa^kLa - H.c] . 
These may be obtained from the time dependent density 
matrix p(t) as: 



(0{t F )) = Tr Op{t F ) 



Tr Oe- lHAtF 



(2) 



For non-equilibrium problems the only approach known 
to be reliable is to compute p(t) by evolving forward from 
an initial condition po as in the second term of Eq. [2] 

We take po = f>Q 0t ® /5 load corresponding to decoupled 
impurity and leads and assume that /5g ot is diagonal in 
the occupation number basis. We evaluate Eq.[2]by writ- 
ing the time evolution operators e A in an interaction 
representation with respect to -f/hyb and expanding pow- 
ers of ifhyb- The bare expansion produces diagrams of 
the form shown in panels (a) and (b) of Fig. [T] The pres- 
ence of two time evolution operators in Eq. [2] means that 
two time contours are required, one running from an ini- 
tial time t = (left side of lower contour) to the measure- 
ment time tp (right side) and the other running back to 
initial time (left side of upper contour, label 2tp indicat- 
ing total time interval along double contour). Hybridiza- 
tion vertices V ka adic ka(7 {V^ aG c\ a<y d,j) occurring at times 
t\...tj arc indicated by heavy (empty) dots as in Ref. [l9| 
and are connected by light lines displaced from the basic 
contour indicating contractions of the lead (c) operators 
computed using p lead and by solid, wavy or dashed lines 
indicating propagation in eigenstates of -ffdot- 

A straightforward evaluation of Eq. [2] thus requires a 
sum over all diagrams, a sum over all contractions of 
lead operators, and an integral over all times. The sign 
problem arising from the powers of i limits diagrammatic 
Monte Carlo studies to situations where the mean pertur- 
bation order is < 10 15- I?} and for this reason only times 
t < 2/r are accessible with this method. We therefore in- 
vestigate a bold method. As an analytic resummation we 
use the non-crossing approximation (NCA) [22], [23| but 
we emphasize that the concepts and methods developed 
here apply to any bold expansion. 

The NCA integrates diagrams without 'crossing' hy- 
bridization lines [Fig. [TJz shows a sample] using coupled 
integral equations. The correction diagrams are obtained 
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FIG. 1. Diagrams arising in the bare hybridization and 'bold' 
expansions of Eq. [TJon the Keldysh contour: propagation in 
eigenstates of H^ot is described by wavy and dashed lines. 
Solid lines denote hybridization functions, and circles hy- 
bridization vertices, a) NCA diagram, b) diagram contain- 
ing crossing, d) bold propagator (dashed double line) which 
resums a), c) Diagram of bold expansion which resums dia- 
grams including b). Diagram e) illustrates vertex resumma- 
tions, and diagram /) shows a bold vertex diagram. 



along lines similar to those of the equilibrium algorithm 
of Ref. [3]. Bare 'atomic state' propagators and non- 
crossing hybridization lines are replaced by 'bold' NCA 
propagators (denoted by heavy lines in Fig.[T]). Diagrams 
with non-crossing hybridization lines contained in the un- 
derlying NCA propagators are not sampled. Thus for ex- 
ample the bold sampling process collapses diagram (a) of 
Fig. Q] to diagram (d), and diagram (b) to diagram (c). 

In defining crossing and non-crossing diagrams the is- 
sue of lines connecting one contour to another arises. 
These lines may be interpreted as vertex corrections to 
the operator placed at the measurement point tp and to 
the initial density matrix. While they can be sampled 
directly, we find that it is best to rcsum non-crossing 
lines spanning a contour into NCA vertex corrections, 
and to replace the bare operator by the combination of 
the operator and its NCA vertex correction. An exam- 
ple of a vertex correction is shown in panel (e) of Fig. [T] 
and the resulting 'bold' version in panel (f). Vertex cor- 
rections, especially for the vertices spanning the initial 
density matrix, significantly reduce the expansion order 
and the dynamic sign problem and allow us to perform 
an expansion about the NCA steady state. 

Our Monte Carlo process is defined by moves which 
propose the addition or removal of vertices on either con- 
tour. The proposals are made without regard to whether 
the diagram is bold or not, but a proposed move which 
produces a diagram which is subsumable into a bold di- 
agram is rejected. The procedure is exact because each 
bare diagram is contained in exactly one bold diagram. 
We combine diagrams in such a way that all terms are 
real 17[ so the phase problem becomes a sign problem. 
For computations of a given observable O the accep- 
tance/rejection probabilities of a given move are deter- 
mined from the absolute value of the contribution to (O) 
and one measures (O sign) / (sign) . The expectation value 
of the sign decays exponentially with perturbation order 
considered and thus with the time interval to be studied. 
Management of the sign problem is a crucial issue in this 
and related methods. 
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FIG. 2. Average sign as a function of time in evaluation of 
expectation value of current for U = 4, /3 = 50, H = 0, V = 5. 
Bare expansion: diamonds. Bold expansion: circles. Other 
lines: Bold expansion truncated at orders 3 to 6. 



Wc have found it useful to define a diagrammatic con- 
figuration at expansion order k (i.e. a set of vertices at 
times t\...t2k) as the sum of all contractions of lead oper- 
ators consistent with the crossing condition. The lack of 
a Wick's theorem means that the sum must be performed 
explicitly. The exponential growth with perturbation or- 
der of the number of possible contractions sets a limit 
~ 10 on the order which can be reached, but this limit 
is less severe in practice than the limit imposed by the 
sign. It is possible that higher orders may be reached 
by integrating diagrams individually or combining only 
a subset of them; this has not yet been explored. 

Fig. [5] shows the time dependence of the expectation 
value of the sign computed in an expansion of the current 
for the noncquilibrium Anderson model. The diamonds 
show the sign obtained from the bare hybridization ex- 
pansion method of Rcfs [l5- 17 1 and the circles the sign 
obtained from a straightforward application of the bold 
method. (Essentially identical sign vs time curves are 
found for all parameters studied except that (sign) in- 
creases at very high T > T.) The larger mean value of 
the sign at a given time in the bold method arises be- 
cause fewer perturbation orders are needed to reach a 
solution. The exponential decrease of (sign) with time 
visible in Fig. [2] constrains the times that can be studied 
with finite resources. We see that the straightforward 
bold expansion can reach w twice as long a time as the 
bare expansion. 

In contrast to equilibrium simulations, where the dia- 
grams generated by the Monte Carlo process are typically 
the ones most important to the evaluation of the observ- 
able, we find that in the noncquilibrium situations con- 
sidered here an unconstrained Monte Carlo exploration 
of bold diagrams generates many high order diagrams 
which sum to zero in the observable. Thus we define a 
Monte Carlo process which considers only diagrams with 
perturbation order less than or equal to some value k, 
and then increase k until convergence is reached. For the 
cases studied here a k < 8 sufficed. Fig. [5] shows that the 
mean sign decreases exponentially with increasing max- 



FIG. 3. Order by order convergence for times t = 1.5, ...,4. 
Upper panel: spin imbalanced parameters, density expansion, 
U = 4, V = 5, state | f) (initial state | blowup to region 
of biggest differences. Lower panel: typical case, U = 8, V = 
2,H — 0, state | "|\|,) starting from | •f). Inset: order-by-order 
contribution to current, for U = 4, V — 5 and times indicated. 



imum perturbation order but for a given perturbation 
order saturates at a non-zero value. We find that once 
the correct k is identified, the bold expansion can be ar- 
ranged so that convergence is uniform in time: the mean 
perturbation order required to obtain a convergent result 
does not increase as the time interval is increased. 

Convergence is poorest for spin-dependent properties 
of spin-imbalanced initial conditions; the main panel of 
Fig-Oshows the slowest-converging case we have encoun- 
tered so far. The straightforward bold expansion only 
converges out to times t sa 2, the convergence with order 
is oscillatory but by 7 th order an acceptable convergence 
is reached as can be seen from the coincidence of the 
6 th , 7 th , and 8 t/l order results. The lower panel shows a 
more typical case, where convergence is monotonic and 
occurs by 4 th order. The upper inset presents the con- 
tribution w(t) made to the current at time t by the sum 
of all diagrams of a given order. We see that for this 
case diagrams of order > 5 make no net contribution, 
but would be extensively sampled in a straightforward 
bold Monte Carlo calculation. 

The much longer times accessible via the methods pro- 
posed here allow us to reach physically interesting steady 
states. Fig. |4] shows the evolution of particular diagonal 
elements of the density matrix for different model param- 
eters and starting from different initial conditions. The 
top and bottom traces (U = 8, V = 1, H = 0.5) show the 
evolution of the spin down (favored by H) and empty 
states starting from the initial condition in which the 
dot is in | J,). The similar time scales in the evolution of 
the empty and singly occupied states show that the time 
dependence, which is rapid and is captured correctly by 
the bare and straightforward bold methods, is almost en- 
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FIG. 4. Time evolution of dot states from specified initial 
condition, for parameters indicated calculated using bare ex- 
pansion, bold expansion, and truncated bold expansion. Dot- 
ted line: NCA. Inset: Decay rate 1/Ti calculated for nonequi- 
librium Anderson model at voltages and temperatures indi- 
cated for {7 = 8. Open symbols: NCA. Filled symbols: Bold. 

tirely due to charge relaxation. By contrast, the middle 
traces (U = 4, V = 5, H = 0) show the evolution of spin 
up and down states from a spin polarized initial condi- 
tion. The much slower spin relaxation is evident. The 
times t > 3 required to access the steady state are only 
accessible by the new methods proposed here. Similarly, 
Fig. [5] shows the time evolution of the current. 

For a wide range of parameters and initial conditions 
we find that the magnetization relaxes exponentially to 
its steady state value: m ~ exp[— t/T\]. The inset of 
Fig. |4] compares the voltage and temperature dependence 
of the spin relaxation rate computed in the bold expan- 
sion and the NCA. The latter systematically underes- 
timates relaxation rates. Remarkably, the temperature 
dependence of T\ is opposite at high and low voltages. 

The presence of analytically computed vertex functions 
in the algorithm allows improved access to the steady 
state by starting the bold computation from the density 
matrix corresponding to the NCA steady state rather 
than from a decoupled or non-interacting initial condi- 
tion. In practice we use the NCA equations to propagate 
forward for a time to from a decoupled state, after which 
the bold interactions are turned on and a further time ts 
is studied. For the parameters we have studied the NCA 
density matrix is typically close to the true steady state 
with the largest differences occurring for a non-zero field 
(dotted line, Fig. 0J. Transients decay quickly. While 
NCA propagators and vertices are required for the en- 
tire time interval tp = to + ts the bold expansion need 
only operate over the much shorter time t b ■ The inset of 
Fig. [5] shows the time evolution of the current from the 
NCA steady state to the numerically exact steady state 
for a representative choice of parameters. The large ini- 
tial transient observed in the main panel is absent. 

In conclusion, we have developed a real time diagram- 



FIG. 5. Time evolution of current for parameters indicated. 
Main panel: starting from empty dot. Inset: starting from 
the NCA steady state. 

matic method that enables a description of long time and 
steady-state properties in nontrivial quantum impurity 
models such as the equilibrium and noncquilibrium An- 
derson models over a wide range of interaction strengths 
and time scales. The approach is based on a systematic 
summation of terms contained in an expansion in powers 
of the hybridization portion of the Hamiltonian about 
a state described by an analytical resummation. This 
'bold expansion' is numerically exact, uniformly conver- 
gent and greatly reduces the real time sign problem that 
inhibits the study of long time properties in 'bare' con- 
tinuous time quantum Monte Carlo methods. We found 
that in many cases the non-crossing approximation pro- 
vides reasonably accurate (< 5%) estimates of the diag- 
onal elements of the steady-state density matrix, but is 
less reliable for relaxation rates. Future work will be de- 
voted to using the method described here as a real-time 
impurity solver for DMFT (where two-time correlators 
are required) and applying it to a range of physically rel- 
evant situations. A crucial question is how far into the 
Kondo regime the method can be pushed. 
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